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Abstract:   The  stability  of  finite  difference  models  of  hyperbolic  partial 
differential  equations  depends  on  how  numerical  waves  propagate  and  reflect 
at  boundaries.   This  paper  presents  an  extended  numerical  example  illustrating 
the  key  points  of  this  theory. 


0.  Iiitrodaction 

In  the  nuincrical  solution  of  hyperbolic  partial  differential  equations  by  finite  differ- 
ences, stability  is  well  known  to  be  a  critical  issue.  As  a  first  step  the  difference  model 
must  satisfy  the  von  Neumann  condition  -  that  is,  the  basic  formula  should  admit  no 
exponentially  growing  Fourier  modes.  For  linear  problems  with  smoothly  varying  coeffi- 
cients and  no  boundaries,  this  is  essentially  the  whole  story,  and  in  fact  if  one  rules  out 
algebraically  as  well  as  exponentially  growing  local  Fourier  modes,  then  stability  is 
assured.  Results  of  this  kind  are  widely  known  and  are  discussed  in  the  superb  book  by 
Richtmyer  and  Morton  [9]. 

When  boundaries  are  introduced,  the  stability  problem  becomes  noorc  subtle.  Even 
here  the  literature  is  copious,  and  a  dozen  or  more  people  have  made  substantial  contribu- 
tions, including  Godunov  and  Ryabenkii,  Strang,  Osher  [7],  Kreiss,  Gustafsson, 
Sundstr5m,  Tadraor,  and  Michelson.  The  best  known  paper  in  this  area  is  the  one  by 
Gustafsson,  Kreiss,  and  SundstrBm  in  1972  [4],  which  presents  what  is  now  often  referred 
to  as  the  "GKS  stability  theory".  The  great  strength  of  the  GKS  paper  h  that  it  estab- 
lishes a  necessary  and  sufficient  stability  condition  for  difference  models  of  very  general 
form  -  three-point  or  multipoint  stencil  in  space,  two-level  or  multilevel  in  time,  explicit 
or  implicit,  dissipative  or  nondissipative.  A  difficulty  with  the  papjer  is  that  it  is  very  hard 
to  read,  and  this  has  regrettably  limited  its  influence.  Fortunately,  some  more  accessible 
accounts  have  appeared  recently,  including  the  report  by  Gustafsson  in  this  volume. 

My  own  work  in  this  field  has  been  concerned  with  giving  the  stability  question  for 
initial  boundary  value  problem  models  a  physical  interpretation  based  on  the  ideas  of 
dispersive  wave  propagation  and  group  velocity.  Group  velocity  effects  in  finite- 
difference  modeling  have  been  surveyed  by  me  in  [9]  and  by  Vichnevetsky  and  Bowles  in 
[13];  others  who  have  been  interested  in  these  matters  include  Matsuno,  Grotjahn,  and 
O'Brien  in  meteorology;  Alfold,  Bamberger,  and  Martincau-Nicoletis  in  geophysics; 
Kentzer,  Giles,  and  Thompkins  in  aerodynamics;  and  Hedstrom  and  Chin  in  theoretical 
numerical  analysis.  I  have  described  a  group  velocity  interpretation  of  the  one-boundary 
stability  problem  in  [10],  showing  in  particular  how  the  GKS  "perturbation  test"  for 
unstable  "generalized  eigensolutions"  is  equivalent  to  a  test  of  the  sign  of  a  group  velo- 
city. In  [11]  this  approach  is  made  rigorous,  and  various  theorems  on  unstable  growth 
rates  are  obtained.  In  [12]  I  have  extended  these  ideas  to  problems  with  two  or  more 
boundaries  or  internal  interfaces,  where  stability  depends  on  what  happens  when  wave 
packets  reflect  back  and  forth.   The  latter  work  was  motivated  in  part  by  ideas  of  Kreiss 
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(see  e.g.  Sec.  7  of  [4]),  of  Beam,  Wanning,  and  Yee  [1],  and  of  Giles  and  Thompkins  [2]. 

There  is  an  analogous  theory  for  p.d.e.'s  rather  than  difference  approximations. 
Again,  wave  radiation  from  the  boundary  is  a  general  mechanism  of  ill-posedness.  See 
Kreiss  [6]  for  the  basic  theory,  and  the  new  survey  by  Higdon  [5]  for  the  wave  interpreta- 
tion. The  difference  is  that  for  p.d.e.'s,  nontrivial  cases  of  ill-posedness  do  not  arise 
unless  the  domain  contains  two  or  more  space  dimensions. 

The  purpose  of  this  paper  is  to  survey  these  results  relating  stability  and  wave  propa- 
gation by  means  of  an  extended  example.  With  the  aid  of  many  illustrations  we  will  see 
exactly  how  waves  can  get  amplified  by  reflection  at  boundaries  and  how  this  can  lead  to 
instability. 
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Figure  1.1.  Dispersion  relation  for  CrEnk-Nicolson  with  X=l. 


Now  that  the  dispersion  function  is  known  we  can  synthesize  v"  as  follows: 
V"  =  ^   f  *'"-f«)'*t«)vO(^)J5,      x=jh,    t  =  nk. 


(1.5) 


Anned  with  this  equation  we  could  duplicate  the  behavior  of  CN  by  computing  Fourier 
integrals.  There  is  little  profit  in  that,  but  what  (1.5)  does  offer  is  the  prospect  of  approx- 
imate evaluation  by  a  stationary  phase  argument.  For  observe  that  the  exponential  term 
introduces  an  oscillatory  behavior  that  will  make  the  integrand  tend  to  cancel  to  zero  if 
<o(^)  and  v^(?)  are  smooth.   The  exception  is  that  at  values  of  ^  satisfying 

there  is  no  oscillation  and  no  cancellation.  In  other  words,  most  of  the  energy  associated 
with  wave  number  ^  travels  approximately  at  the  group  velocity 
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C  =  - 


dk 


(1.6) 


Of  course  this  argument  is  vague,  but  it  can  be  made  precise  in  various  ways;  see  for 
example  Lemma  5.1  of  [11]. 

Thus  wave  energy  travels  at  a  velocity  given  by  the  negative  of  the  slope  of  the 
dispersion  relation.   For  CN  we  can  differentiate  (1.4)  implicitly  to  obtain 


-AS  .3  1  ,•(1-)   :. 
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C  =  -cos^  cos 


20)* 


(1.7) 


EUminating  w  by  means  of  (1.4'),  or  differentiating  (1.4')  directly,  converts  this  to  the 
functional  fonn 


C(4)  = 


-cos£,h 


1  +  -^sin=5A 
4 


(l.T) 


This  function  is  plotted  in  Fig.  1.2.  One  sees  that  for  well-resolved  waves  -  i.e.  ^«0, 
or  many  points  per  wavelength  -  energy  travels  at  velocity  -1,  as  it  should  according  to 
the  p.d.e.  u,  =  «,.  Less  well  resolved  waves  have  lower  speeds  (less  negative  velocities), 
and  it  is  this  fact  that  gives  rise  to  familiar  oscillations  around  discontinuities.  At  the 
extreme,  the  sawtoothed  (or  parasitic)  wave  vj  =  (-ly,  i.e.  ^  =  ±'n,  has  group  velocity 
+  1,  so  energy  in  this  mode  travels  in  the  physically  wrong  direction. 
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Figure  1.2.   Group  velocity. 


Let  us  confirm  this  last  prwliction  by  an  experiment.    Fig.  U  shows  the  evolution 
under  CN  with  h  =  .01  and  X  =  1  of  an  initial  signal 

v;  =  [i+(-iy]  *-*"<•- 5)^ 

on  the  interval  [0,1].    This  "rectified  Gaussian",  shown  in  Fig.   1.3a,  contains  equal 
amounts  of  energy  at  |/i  =  0  and  at  ^  =  ±tt.   Figures  1.3b  and  1  Jc  show  the  wave  forms 
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at  times  .2  and  .4.  As  predicted,  the  two  wave  (ximponcnts  have  separated  and  traveled 
in  opposite  directions.  This  backwards  motion  of  the  parasitic  wave  component  is  of 
course  a  purely  numerical  effect.   For  further  examples  sec  [9]  and  [13]. 


(a) 


f=0 


(b) 


r=.2 


(c) 


i\Vv    t=A 


Figure  1.3.     Propagation  of  energy  at  group  velocities  C=  ±1. 


In  analyzing  the  behavior  of  an  arbitrary  difference  model,  a  useful  question  to  ask 
is:  given  a  frequency  (Dq,  what  associated  wave  numbers  4  are  admitted  by  the  dispersion 
relation,  and  how  do  the  corresponding  waves  exp(i(a)o/+^))  propagate?  To  get  the 
answer  for  CN,  imagine  drawing  a  horizontal  line  at  height  coq  in  Fig.  1.1.  Assuming  wq 
is  small  enough,  it  will  intersect  the  dispersion  curve  at  two  values  ^j  and  ^j-  Of  the  two 
corresponding  sine  waves,  one  has  C^O  and  one  has  CaO.  For  simplicity,  from  now  on 
we  will  write  t|  for  the  wave  number  corresponding  to  C^O  and  ^  for  the  other  one.  By 
(1.4),  5  and  r\  are  related  under  CN  by 


?:  rightgoing 
•q:  leftgoing 


(1.8) 


The  same  relationship  holds  for  any  difference  model,  such  as  leap  frog  or  backwards 
Euler,  whose  spatial  discretization  consists  of  the  usual  second-orxfer  centered  difference. 
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As  a  more  complicated  example,  suppose  we  had  a  nondissipative  difference  formula 
with  the  dispersion  function  plotted  in  Fig.  1.4.  (An  arbitrary  continuous  dispersion  func- 
tion  defines  a  bounded  operator  vVv"*^  in  /  ,  a  Fourier  multiplier,  but  this  operator  can  be 
realized  by  finite  differences  only  when  the  function  is  the  solution  of  a  trigonometric 
polynomial  equation  in  u  and  5.)  According  to  the  figure,  there  arc  four  wave  numbers  ^ 
associated  with  the  frequency  (Oj.  gj  and  ^3  correspond  to  leftgoing  waves,  and  ^2  *^d  tt 
to  rightgoing  waves. 
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Figure  1.4.   Dispersion  relation  for  an  unknown  difference  formula. 


On  the  face  of  it  the  situation  looks  different  at  frequency  0)2  ~  there  is  only  one  wave 
propagating  in  each  direction.  However  in  fact,  the  missing  two  wave  numbers  still  exist, 
but  they  have  become  complex.  One  has  negative  imaginary  part,  and  corresjxmds  to  a 
leftgoing  evanescent  wave  that  does  not  carry  energy;  the  other  has  positive  imaginary 
part  and  is  evanescent  and  rightgoing. 

Actually,  as  far  as  one  can  tell  from  Fig.  1.4,  evanescent  modes  may  exist  at  fre- 
quency (Ji  too.  TTiis  depends  on  the  size  of  the  stencil  of  the  difference  formula.  The 
general  situation  is  this:  for  a  scalar  difference  formula  with  stencil  extending  /  points  to 
the  left  of  center  and  r  points  to  the  right,  the  dispersion  relation  is  a  trigonometric  poly- 
nomial in  ^  of  order  l+r,  which  for  each  u  has  l+r  complex  roots  {  that  break  down  into 
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exactly  r  "leftgoing"  and  /  "rightgoing"  linearly  independent  wave  modes.   Here  we  say 
that  a  wave  exp(j(ujr+4r))  with  lmu)  =  0  (or  more  generally  ImwsO)  is  rightgoing  if  either 

(l)Im4  =  0andC(?,a))2:0 

or 

(2)  lm^>0. 

A  Uftgoing  wave  is  defined  with  the  directions  of  the  inequalities  in  (1)  and  (2)  reversed. 
These  definitions  and  their  consequences  are  studied  in  detail  in  [11]. 

Thus  the  general  solution  of  the  form  yj"  =  f'"'<J>y  to  an  arbitrary  scalar  finite  differ- 
ence formula  can  be  written 


v;   =   *"-'2ay^'  +  ^'"'2:^v'"^.        x=ih,    r=nJt. 

v-l  v-1 

(rightgoing)            (leftgoing)  \^-^f 

For  convenience  I  have  relabeled  a^+i a;+^  by  pj p^,  and  g^+j ^,^^  by 

"Hi.   •  •  •  .^r' 


-:ivbnt?oa  nfj  i\;r-A 


Ti      ^^obrb' 


"C-- 
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2.  One  boundary:  reflectioo  coefficients 

What  if  we  new  let  more  time  elapse  in  Fig.  1.3,  so  that  the  waves  hit  the  boun- 
daries at  x=  ±1?  The  result  will  of  course  depend  on  the  boundary  conditions  there. 

The  general  principle  for  analyzing  such  problems  is  to  look  for  solutions  containing 
only  a  single  frequency  coq.  Different  frequencies  can  be  superposed  later.  If  a  wave 
txp(i(<t>(^+^))  hits  a  boundary,  then  the  reflected  result  after  the  initial  transient  has  died 
away  will  be  a  linear  combination  of  waves  ^a.^,exp(i(ii(jt+^^))  with  the  same  frequency  uq 
but  with  various  new  wave  numbers  4^.  In  general  the  wave  numbers  in  the  reflected 
waves  will  be  all  of  those  fulfilling  the  following  two  conditions: 

(1)  wq.Cv  satisfy  the  dispersion  relation. 

(2)  The  wave  exp(i(<iiQi+i^^))  propagates  away  from  the  boundary  into  the  interior  (the 
radiation  condition).  This  means  that  any  wave  reflected  at  a  lefthand  boundary 
must  be  rightgoing,  and  any  wave  reflected  at  a  righthand  boundary  must  be  leftgo- 
ing. 

Neither  of  these  statements  mentions  the  boundary  conditions.  Those  do  not  affect  the 
set  of  reflected  waves,  just  the  coefficients  a^,. 

To  compute  numerical  reflection  coefficients  one  simply  inserts  the  wave  (1.9)  into 
the  numerical  boundary  conditions.  For  a  general  formulation  consider  a  lefthand  boun- 
dary at  x=j=0,  and  let  the  numerical  boundary  conditions  there  consist  of  /  linear  homo- 
geneous equations  giving  v^*^ vj**/  as  linear  combinations  of  other  values  v7.  Inser- 
tion of  (1.9)  yields  a  linear  reflection  equation 

£(u))Qt  =  D(u)3  (2.1) 

relating  leftgoing  and  rightgoing  wave  coefficients.  Here  E  and  D  are  matrices  of  dimen- 
sions /x/  and  /xr,  respectively.  If  £(a))  is  nonsingular  we  can  solve  for  the  rightgoing 
wave  coefficients  to  get 

a  =  ^(a))3  =  [£(u))]-^D(a>)3.  (2.2) 

A(oi)  is  called  the  reflection  matrix.  Compare  [4],  cq.  (10.2). 

Figure  2.1  shows  what  happens  when  the  integration  of  Fig.  1.3  is  carried  to  r>iy2 
with  the  following  conditions  at  the  boundaries: 

v2*i  =  v7^i,  (2.3a) 

v7*i  =  0.  (2.4a) 


1  ■     r^ri't^i^ -rV""r    nf"" 
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Herc  J=l/h  is  the  index  of  the  grid  point  at  the  boundary  x=l.  Figs.  2.1a  and  2.1b  show 
the  configuration  at  time  r=.6  and  r=.8. 


r=.6 


/=.8 


Figure  2.1.   Continuation  of  Fig.  1.3  showing  reflection  at  boundaries. 


Clearly  the  leftgoing  pulse  with  ti=0  has  generated  a  rightgoing  reflected  wave  with 
^=v/h  at  very  small  amplitude.  The  rightgoing  pulse  with  ^=v/h  has  generated  a  consid- 
erably larger  reflection  with  t|  =  0.   Let  us  predict  these  reflected  amplitudes  after  the  fact. 

First  we  compute  the  reflection  coefficient  for  (2.3a).    In  the  present  case  (1.8) 
implies  that  (1.9)  reduces  to 

v;  =  [a<r't'  +  3<r'n']«''"  =  [(-iya*-"''  +  p*""]<f'"',       x=jh,   t=nk.  (2.5) 

For  simplicity  in  such  computations  it  is  convenient  to  introduce  the  abbreviations 
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K  =  e'^,      p.  =  e'^,     z  =  e" 
Then  (2.5)  becomes 

vy"  =  (aK>'  +  (3n/)z". 
For  CN  (1.8)  becomes  K  =  -l/^I.   Inserting  (2.7)  in  (2.3a)  gives 

a+P  =  (oK  +  Pp.)  =  -o/ji+Pti. 
In  the  form  of  (2.1)  this  is 


(2.6) 


(2.7) 


that  is, 


(l  +  ^)a  =  (ti-l)p. 


3       %+l  2 


(2Jb) 


(23c) 


.      »-^^jL/i 


.;^lt.  cni02  •\d)!!nic.. 
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For  -11=0,  as  in  the  present  experiment,  we  get  a/p~0,  and  this  explains  the  very  small 
amplitude  of  the  pulse  reflected  from  the  lefthand  boundary  in  Rg.  2.1. 

Now  the  analogous  computation  for  (2.4a).   Inserting  (2.7)  gives 

K-'a  =  -pL-'P,  (2.4b) 

that  is, 

|-  =  -(-jiY  =  i-\y^V-^\  (2.4c) 

Thus  the  amplitude  of  the  reflected  wave  should  be  equal  and  opposite  to  that  of  the 
incident  wave,  regardless  of  tj.  Fig.  2.1  confirms  this  nicely. 

Let  us  return  to  the  lefthand  boundary  and  consider  some  alternative  boundary  con- 
ditions.  Suppose  we  replace  (2.3a)  by 

vr^  =  v^.  (2.8a) 

Then  insertion  of  (2.7)  leads  to 

iz-^-)a  =  i-z+y.)^,  (2.8b) 

or 

f  =  y^f^  =  .•  .'V-sin^l^/cos^:!^.  (2.8c) 

Again  this  predicts  a  near-zero  reflected  amplitude  forT|  =  0~(D.   On  the  other  hand  sup- 
pose we  impose 

v^^  =  v^^.  (2.9a) 

The  reflection  equation  is  then 

(1-1/ji^a  =  (^2-1)3.  (2.9b) 

This  implies 

I"  =  ti-  =  f^^  (2.9c) 

much  as  in  (2.4c),  at  least  for  ii^stO.   At  ti=0  (2.9b)  has  the  form  Oa=Op,  so  it  is  not  solv- 
able except  in  a  limiting  sense. 

Figure  2.2  plots  results  of  experiments  with  boundary  conditions  (2.8a)  and  (2.9a). 
Fig  2.2a  shows  the  initial  condition,  a  leftgoing  wave  with  ti=0  on  a  mesh  with  A=.01  for 
[0,1],  and  Figs.   2.2b  and  2.2c  show  the  results  at  /=1   under  (2.8a)  and  (2.9a), 
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respectively.  The  predicted  reflection  coefficients  0  and  1  are  clearly  in  evidence.  In  fact 
no  reflected  energy  at  all  is  visible  in  Fig.  2.2a;  this  is  because  with  X=l,  i.e.  h=k,  (2.8c) 
is  identically  zero. 


20 
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f{v^ 


(2.9a) 


r=l 
(2.10a) 


Figure  2.2.   Effect  of  various  lefthand  boundary  conditions. 


Figure  2. 2d  shows  the  response  of  CN  to  the  same  initial  data  with  a  more  contrived 
boundary  condition: 

v2*i  =  -v;-*!  +  vf  1  +  vf  1.  (2.10a) 

This  time  a  marked  instability  is  evident  (note  the  vertical  scale).    To  see  why,  compute 
the  reflection  equation 

which  can  be  simplified  to 

(l+»i)(l-»x)2a  =  -^t'(l  +  n)2(l-»jL)p.  (2.10b) 

that  is, 


i  =   -a'l±ii   =  i^'r^COX^. 


3 


1-n 


(2.10c) 


■";-?:: 
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At  71=0,  this  reflection  coefficient  is  infinite. 

The  existence  of  a  frequency  w  at  which  the  reflection  coefficient  at  a  boundary  is 
infinite  implies  that  a  finite-difference  model  is  unstable.  In  [11]  it  is  shown  that  under 
reasonable  assumptions,  such  a  model  potentially  amplifies  initial  data  in  the  /*  norm  at  a 
rate  at  least  proportional  to  the  time  step  number  n.  Obviously  this  is  an  unstable  situa- 
tion, since  as  h  and  k  are  decreased,  the  index  n  corresponding  to  a  fixed  time  t  increases 
to  ». 

The  presence  of  an  infinite  reflection  coefficient  is  a  stronger  condition  than  "GKS- 
instability,"  a  name  sometimes  given  to  instability  according  to  Defn.  3.3  of  [4].  In  fact  a 
difference  nxxfel  is  GKS-unstable  if  and  only  if  the  matrix  £((i))  in  (2.1)  is  singular  for 
some  (1)  with  Imu^O,  i.e.  |z|al  (Thm.  la  of  [11]).  For  scalar  problems  with  three-point 
stencils  this  amounts  to  the  condition  that  the  coefficient  of  a  on  the  left  side  of  a  reflec- 
tion equation  such  as  (2.3b),  (2.4b),  (2.8b),  (2.9b),  or  (2.10b)  is  0  for  some  w.  Thus,  for 
example,  the  boundary  condition  (2.9a)  is  GKS-unstable  for  CN,  since  the  Icfthand  side 
of  (2.9b)  is  zero  at  ji,=z=l,  even  though  no  apparent  catastrophe  occurred  in  Fig.  2.2c. 
In  such  situations  the  reflection  coefficient  remains  finite  if  it  happens  that  the  righthand 
side  of  the  reflection  equation  has  a  zero  at  the  s  me  value  of  a>  of  at  least  as  high  an 
order.  When  this  occurs,  it  is  shown  in  [11]  that  unstable  growth  of  initial  data  in  the  /*■ 
norm  need  proceed  no  faster  than  in  proportion  to  Vn. 

To  summarize,  there  are  at  least  three  distinct  drcumstanoes  that  may  obtain  at  the 
boundary: 

(1)  nonsingular  refl.  eq.,  finite  refl.  coeff.:  stable  (  ||v||=0(l)  ); 

(2)  singular  refl.  eq.,  finite  refl.  coeff.:  weakly  unstable  (  ||v||=0(V^)  ); 

(3)  singular  refl.  eq.,  infinite  refl.  coeff.:  unstable  (  ||v||=0(n)  ). 
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3.  Two  boundaries:  reflection  back  and  forth 

What  if  wc  let  even  more  time  elapse  in  Figs.  1.3  and  2.1?  The  answer  is  simple: 
new  reflections  will  occur  as  the  two  wave  packets  bounce  back  and  forth  between  ^=0 
and  x=l.  Each  time  a  rightgoing  wave  packet  with  i~v/h  hits  x=l,  it  will  reflect  as  a 
leftgoing  wave  packet  with  11=0  having  essentially  the  same  amplitude.  But  each  time  a 
leftgoing  wave  hits  j:=0,  it  will  reflect  as  a  rightgoing  wave  of  greatly  diminished  ampli- 
tude. The  total  energy  in  [0,1]  will  consequently  decay  rapidly  to  near  0.  The  reason  one 
can  only  say  "near"  is  that  a  certain  portion  of  the  energy  will  have  5~tt/2A,  and  since 
C=0  at  this  wave  number,  this  portion  tends  to  stay  in  a  fixed  position  and  never  hit  the 
boundaries. 

Figs.  3.1a  and  3.1b  confirm  this  prediction  by  showing  the  numerical  solutions  at 
f=2  and  f=4.  At  r=2  the  signal  that  remains  is  very  weak,  and  at  r=4  it  is  somewhat 
weaker. 


r=2 


t=A 


Figure  3.1.   Further  continuation  of  Figs.  1.3  and  2.1. 


Suppose,  more  generally,  we  set  up  the  CN  model  on  an  interval  [0^]  with  some  pair 
of  stable  boundary  conditions  at  x=0  and  x=L  for  which  the  reflection  coefficient  func- 
tions are  A^(ui)  and  ^2(0)),  respectively.  If  v°  is  a  wave  packet  with  frequency  u,  then  as  t 
increases  this  packet  will  reflect  back  and  forth,  undergoing  amplification  by  [Aj  or  IA2I 
each  time  it  hits  a  boundary.  Let  C(a))  denote  the  group  speed  |C(5,(d)|  =  |C(Ti,a))|.  Then 
the  travel  time  for  a  complete  circuit  involving  a  reflection  at  each  boundary  will  be 

2L 


n«-)  = 


C(a.) 


(3.1) 
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After  such  a  circuit  the  amplitude  will  have  increased  by  a  factor 

a(a))  =  lAi(u,>42M|.  (3.2) 

Combining  these  formulas,  we  reach  the  following  conclusion:  as  t-»,  energy  associated 
with  frequency  o)  grows  approximately  at  the  rate 

Jj^  -  a(u,)^^«-)  =  aCa,)^^^)^.  (3.3) 

For  o>l  this  equation  predicts  exponential  growth,  and  for  a<l,  exponential  decay. 

On  the  basis  of  this  analysis  one  can  now  make  some  very  interesting  observations 
about  stability  and  convergence.  The  following  ideas  are  spcUed  out  in  detail  in  [12]. 

The  first  observation  is,  stability  does  not  preclude  exponential  growth  in  time.  Wc 
have  seen  in  the  last  section  that  stability  for  a  single  boundary  implies  A<<»,  but  not 
i4sl.  Thus  it  is  easily  possible  for  a  model  composed  of  two  individually  stable  boun- 
daries to  generate  growth  of  truncation  errors  at  a  rate  cxp(constr).  One  might  think  that 
this  reveals  that  the  concatenation  of  two  stable  boundaries  is  in  general  unstable.  But  in 
fact,  exponential  growth  does  not  imply  instability  in  the  usual  Lax-Richtmyer  sense,  and 
it  will  not  prevent  convergence.  The  reason  is  that  the  growth  factor  for  any  fixed  t  does 
not  get  larger  as  h,k-^,  whereas  the  truncation  errors  that  the  factor  multiplies  get  smaller 
(by  consistency). 

These  conclusions  arc  already  known.  In  fact  Thm.  5.4  of  [4]  asserts  that  in  general, 
the  concatenation  of  two  GKS-stable  difference  schemes  is  always  GKS-stable,  while  Sec 
7  of  [4]  investigates  the  conditions  under  which  exponential  growth  in  t  will  occur  for  a 
particular  stable  2x2  example.  What  is  new  here  is  the  interpretation  by  reflection  coeffi- 
cients. This  interpretation  also  helps  to  clear  up  a  misconception.  Some  people  have 
thought  that  although  exponential  growth  can  occur  for  finite  h  and  k,  it  must  vanish  as 
h,k-0.  But  the  study  of  reflecting  wave  packets  gives  no  reason  to  expect  such  good  for- 
tune, and  an  example  in  [12]  proves  that  indeed  it  Is  not  to  be  expected.  To  be  sure  of 
eliminating  growth  by  refining  the  mesh,  one  needs  a  difference  formula  that  is  suitably 
dissipative  [3]. 

There  is  a  reason  why  expxmential  growth  in  t  may  be  a  problem:  in  practice,  time 
dependent  finite  difference  models  arc  often  used  (especially  in  aerodynamics)  to  compute 
approximate  solutions  for  the  steady  state  r-«.  Obviously  growth  of  errors  at  a  rate  e°*^ 
will  make  such  a  procedure  fail.  With  this  in  mind.  Beam,  Warming,  and  Yec  [1]  have 
defined  a  difference  model  to  be  P -stable  if  it  is  GKS-stabk,  and  in  addition,  it  admits  no 
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exponentially  growing  modes.   In  this  language  the  above  observation  becomes:  stability 
does  not  imply  P-stability. 

Our  second  observation  now  comes  as  the  natural  complement  to  the  first:  having 
reflection  coefficients  bounded  by  1  does  preclude  exponential  growth  in  time.  For  in  (3.3) 
we  see  that  if  a(a))<l  for  every  u,  then  no  frequency  exists  that  can  experience  repeated 
amplification  by  reflection.  This  argument  can  be  made  rigorous,  and  Prop.  11  of  [12] 
states:  if  a(w)<l  for  every  o),  a  difference  model  is  P-stablc. 

For  a  third  and  final  observation,  let  us  turn  to  the  case  in  which  one  or  both  boun- 
daries is  individually  unstable.  Then  we  find,  a  mildly  unstable  boundary  with  an  infinite 
reflection  coefficient  may  become  catastrophically  unstable  when  a  second  boundary  is  intro- 
duced. For  consider  a  boundary  condition  which  has  A(cx>q)=»  for  some  uq.  If  a  wave 
packet  at  frequency  wq  hits  this  boundary,  it  will  be  amplified  by  a  large  factor  (typically 
0{S),  where  //  is  the  width  of  the  packet  in  grid  points).  For  a  single  boundary,  this  is  a 
one-time-only  event,  as  in  Fig.  2. 2d.  But  when  there  are  two  boundaries,  the  reflected 
wave  may  reflect  at  the  o±er  boundary,  then  return  to  be  amplified  again,  and  so  on. 
For  boundaries  separated  by  N  grid  points  one  gets  growth  potentially  at  a  very  rapid  rate 

||v(0)||        ^  "  • 

This  instability  is  a  far  cry  from  the  growth  proportional  to  n  mentioned  in  the  last  sec- 
tion, and  it  renders  the  difference  model  totally  useless. 

This  phenomenon  of  catastrophic  two-boundary  interactions  was  first  noted  long  ago 
by  Kreiss.  The  advantage  of  the  present  point  of  view  is  that  it  makes  it  clear  thst  the 
problem  is  associated  not  with  all  unstable  boundaries,  but  only  with  those  having  infinite 
reflection  coefficients.  Unstable  boundary  conditions  with  finite  reflection  coefficients 
typically  exhibit  only  weak  instability  even  when  other  boundaries  are  present.  Of  course, 
sometimes  weak  instabilities  may  be  more  dangerous  than  strong  ones,  since  they  can 
more  easily  go  undetected. 

Our  final  numerical  experiment,  sununarizcd  in  Fig.  3.2,  illustrates  the  difference  in 
two-boundary  behavior  between  unstable  boundary  conditions  with  finite  and  infinite 
reflection  coefficients.  Fig.  3.2a  shows  the  initial  data,  a  uniformly  distributed  random 
signal  on  the  usual  CN  mesh  on  [0,1].  At  x=l  the  boundary  condition  is  (2.4a).  First  an 
integration  was  performed  with  the  boundary  condition  (2.9a)  at  x=0,  which  is  unstable 
but  by  (2.9c)  has  a  finite  reflection  coefficient.  In  this  process  nothing  dramatic  occurs, 
no  matter  how  large  t  is.    Fig.  3.2b  shows  the  result  at  f=10,  and  it  looks  qualitatively 
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(2.10a) 


Figure  3.2.    Two-boundary  interactions 


much  like  the  initial  data.  Of  course  this  does  not  imply  that  this  model  will  give  accurate 
answers  to  physical  problems,  but  it  does  reveal  that  GKS-instability  of  a  boundary  condi- 
tion does  not  by  itself  lead  to  explosive  two-boundary  interactions. 

In  contrast,  Rgs.  3.2c-d  show  results  obtained  with  boundary  condition  (2.10a), 
which  has  an  infinite  reflection  coefficient.  The  plots  show  times  r=l,3.  Note  the  vertical 
scales. 


In  summary,  for  two- boundary  problems  one  has  the  following  possibilities: 

(1)  stable  b.c's,  Osa^l:  stable  and  P-subIc  (|lvl|=0(l)) 

(2)  stable  b.c's,  l<o<«:  stable  but  P-unstable  (|Iv||=C?(if'^")) 

(3)  unstable  b.c's,  Oso<ac:  weakly  unstable  (||v||=0(V^)) 

(4)  unstable  b.c's,  a=":  catastrophically  unstable  (}\v\\=0(N°^^)) 
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Before  closiag,  a  word  must  be  said  about  dissipation.  I  have  presented  a  picture  in 
which  numerical  waves  bounce  back  and  forth  between  boundaries,  with  changes  of 
amplitude  occurring  at  each  reflection  but  not  during  the  transit  in  between.  Under  a  dis- 
sipativc  difference  formula,  however,  any  wave  with  5^0  decays  as  it  travels.  The  result 
is  that  the  two-boundary  interactions  I  have  spoken  of  rarely  take  place,  so  that  in  prac- 
tice, the  behavior  of  a  dissipative  two-boundary  model  is  not  much  different  from  what 
the  individual  boundaries  would  suggest.  The  exceptions  occur  when  the  dissipation  is 
very  weak,  as  in  the  problems  with  large  mesh  ratios  considered  in  [1];  or  when  the 
number  of  grid  points  between  boundaries  is  very  small,  as  in  certain  cases  discussed  in 
Sec.  2  of  [12]. 

For  an  overall  summary,  one  may  say  that  difference  models  for  hyperbolic  partial 
differential  equations  exhibit  three  important  physical  processes:  propagation,  reflection, 
and  dissipation  of  waves.  In  general  the  waves  involved  are  numerical  objects  with  no 
physical  meaning,  but  they  are  still  important  to  stability.  An  exact  mathematical  analysis 
of  a  finite  difference  model  is  usually  too  difficult  to  be  practical,  but  even  for  compli- 
cated models  one  can  often  gain  insight  fairly  easily  by  considering  these  physical 
processes. 
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